%% SCRIPT TO EXTRACT SYNTHETIC RADAR WIND OBSERVATIONS FROM WRF GRIDDED DATA DURING THE MARSHALL FIRE

close all
clear all
addpath('~/Dropbox/MATLAB')
addpath('~/Dropbox/MATLAB/SGP_LIDAR_CODES/');

%% Load Marshall radar plume data
load('marhsall_plume_d1.mat');
xv=unique(xxx(:)); %construct an x vector
zv=unique(zzz(:)); %construct a y vector
alt=1; %double check this later...
lons=stationlon; % change variable name so it doesn't get overwritten
lats=stationlat; %change variable name so it doesn't get overwritten
re=6.37*10^6; %earth radius
scalefactor=(re.*cosd(nanmean(lats))).*(pi./180); %scaling for lons<->x distances

%% LOAD WRF DATA FILES
fname='/Users/nlareau/Dropbox/MARSHALL_FIRE/wrfout_d02_2021-12-30_22:15:00';
ncdisp(fname)
wtime=ncread(fname,'Times'); %read character string of the simulation time
T=ncread(fname,'T'); %read the temperature data
U=double(ncread(fname,'U')); %u winds
Uus=(U(2:1:end,:,:)+U(1:1:end-1,:,:))./2; %destagger u grid

V=double(ncread(fname,'V')); %v winds
Vus=(V(:,2:1:end,:)+V(:,1:1:end-1,:))./2; %destagger v grid
W=double(ncread(fname,'W')); %w winds
Wus=(W(:,:,2:1:end)+W(:,:,1:1:end-1))./2; %destagger w grid

%% examine surface map
latsW=squeeze(double(ncread(fname,'XLAT'))); %lat (y) grid for the w data points
lonsW=squeeze(double(ncread(fname,'XLONG'))); %lon (x) grid for the w data points

thgt=squeeze(double(ncread(fname,'HGT'))); %terrain height field
figure(1);clf;
subplot(1,2,1);cla; hold on;
pcolor(lonsW,latsW,squeeze(U(1:end-1,:,1)));shading flat; caxis([-10 20])


hold on;

plot(stationlon,stationlat,'*k');
plot(lonsW(:,300),latsW(:,300),'--k')
colormap(rbmapper(20,-10))
daspect([1 1 1])
caxis([-10 20])
subplot(1,2,2);cla; hold on;
pcolor(lonsW,latsW,squeeze(V(:,1:end-1,1)));shading flat;
hold on;
plot(stationlon,stationlat,'*k');
plot(lonsW(:,300),latsW(:,300),'--k')
colormap(rbmapper(10,-10))
daspect([1 1 1])
caxis([-10 10])

%% compute the meridional variance of the u wind across the marshall fire domain
usfc=Uus(:,:,1);
vsfc=Vus(:,:,1);

uvar=nanstd(usfc(:,100:350),[],2);
vvar=nanstd(vsfc(:,100:350),[],2);

ubar=nanmean(usfc(:,100:350),2);
figure(10);clf;hold on;set(gcf,'color','w')
sp1=subplot(3,1,1:2);cla; hold on;
pcolor(lonsW,latsW,usfc);shading flat;caxis([-10 30])
colormap(rbmapper(30,-10));
plot(lonsW(:,350),latsW(:,350),'--k','linewidth',2)
plot(lonsW(:,150),latsW(:,150),'--k','linewidth',2)
plot(stationlon,stationlat,'*k','markersize',10,'linewidth',2);

daspect([1 .78 1]);
hold on;
xlim([-105.45 -105])
ylim([39.79 40.1])
box on
cbh=colorbar;ylabel(cbh,'m/s')
set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')

%% compute local variance
%across wind std
f1=ones(1,25);
uxstd=stdfilt(usfc,f1);
%along wind std
f2=ones(25,1);
uastd=stdfilt(usfc,f2);


%% add in the marhsall fire burn scar
fileName='marshall.json';
str = fileread(fileName); % dedicated for reading files as text
test=jsondecode(str);
for aa=1:size(test.features.geometry.coordinates)

    pnow=[squeeze(test.features.geometry.coordinates{aa})];
    plot(sp1,pnow(:,1),pnow(:,2),'m','linewidth',2)
end


subplot(3,1,3);cla; hold on;
plot(lonsW(:,1),((uvar)),'r','linewidth',2)
plot(lonsW(:,1),((vvar)),'b','linewidth',2)
plot(lonsW(:,1),ubar,'-k','linewidth',2)
xlim([-105.45 -105])
box on
cbh=colorbar;cbh.Visible='off';
ylabel('m/s')
legend('std(u)','std(v)','Ubar')
set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')



%% Contruct geopotential data
gph=(squeeze(double(ncread(fname,'PH')))+squeeze(double(ncread(fname,'PHB'))))./9.81;
gphus=(gph(:,:,2:1:end)+gph(:,:,1:1:end-1))./2; %destagger geopotential grid

%% Now compute radial components of the wind relative to the radar site

%start by constructing a matrix that contains the angle between a given WRF
%point and the radar location
azimuths=azimuth(stationlat,stationlon,latsW,lonsW);
distancer=sqrt(((stationlon-lonsW).*scalefactor).^2 + ((stationlat-latsW).*111180).^2);
edist=atand((gph-stationalt)./distancer);

figure(2);clf;
subplot(1,2,1);cla; hold on;
% pcolor(lonsW,latsW,azimuths);shading flat;
%   pcolor(lonsW,latsW,distancer);shading flat;
eds=squeeze(edist(:,:,1));
pcolor(lonsW,latsW,squeeze(edist(:,:,1)));shading flat;
plot(stationlon,stationlat,'*k');

caxis([-4 4])

subplot(1,2,2);cla; hold on;
vr=(squeeze(U(1:end-1,:,1)).*sind(azimuths))+ (squeeze(V(:,1:end-1,1)).*cosd(azimuths));
vrcopy=vr; vrcopy(distancer>1*10^4)=nan
pcolor(lonsW,latsW,vrcopy);shading flat;
colormap(rbmapper(10,-10))
caxis([-10 10])
plot(stationlon,stationlat,'*k');
plot(lonsW(:,300),latsW(:,300),'--k')

%% NEXT DO IT FOR REAL

%step one create an interpolant for U wind
%convert lat/long to x-y distances relative to radar
xlon=(stationlon-lonsW).*scalefactor; %lons-> xdistances from the radar
ylat=(stationlat-latsW).*111180; %lats-> ydistances from the radar
xlm=repmat(xlon,1,1,size(U,3)); %create a 3D matrix of these values
ylm=repmat(ylat,1,1,size(U,3)); %create a 3D matrix of these values
dmat=repmat(distancer,1,1,size(U,3));

Utrim=squeeze(Uus);
Vtrim=squeeze(Vus);
Ztrim=squeeze(gphus);
FU=TriScatteredInterp(xlm(dmat<3*10^4),ylm(dmat<3*10^4),Ztrim(dmat<3*10^4),Utrim(dmat<3*10^4));
FV=TriScatteredInterp(xlm(dmat<3*10^4),ylm(dmat<3*10^4),Ztrim(dmat<3*10^4),Vtrim(dmat<3*10^4));

%% Now read in a scan from the DOW and interpolate U to the DOW polar grid points

files=dir('./deployment1/cfrad*.nc');

for ff=288%1:length(files);%50:170%140:170%158;%159:170%1:length(files)

    fname=strcat(files(ff).folder,'/',files(ff).name);
    % fname='cfrad.20211230_222630.637_to_20211230_222642.482_DOW6low_SUR.nc'
    %     ncdisp(fname)

    lats=ncread(fname,'latitude');
    lons=ncread(fname,'longitude');
    alt=ncread(fname,'altitude');
    vel=double(ncread(fname,'VEL'));
    velF=double(ncread(fname,'VEL_F'));
    dbzhc=double(ncread(fname,'DBZHC'));
    azimuth=double(ncread(fname,'azimuth'));
    elevation=double(ncread(fname,'elevation'));
    ranger=double(ncread(fname,'range'));
    re=6.37*10^6;
    ts=double(ncread(fname,'time'));%time since volumne start
    starttime=char(ncread(fname,'time_coverage_start'));
    starttime=[starttime'];
    starttime(1,11)=' ';starttime(20:end)='';
    timer(ff)=datenum(starttime,'yyyy-mm-dd HH:MM:SS');

    clear xdist ydist dbznow

    for ev=1%:size(elevation)

        %%
        evnow=mean(elevation); % pull out the elevation for this scan.
        if evnow>2 && evnow<4.3
            ae=(4/3).*(6.37*10^6); %earth radius
            z=(((ranger.^2) + (ae.^2) + (2*ranger.*ae.*(sind(evnow)))).^(1/2))-ae;
            hdist=  ae.*(asin((ranger.*cosd(evnow))./(ae+z)));%ae.*asin(rdist_HI.*(cosd(evnow))./(ae+z));
            [azmt,hdstmt]=meshgrid(azimuth,hdist); % create a mesh of azimuths and along ground distances
            [~,zdist]=meshgrid(azimuth,z);
            xdist=sind(azmt).*hdstmt; %x component of the horizontal distance
            ydist=cosd(azmt).*hdstmt;
            scalefactor=(re.*cosd(nanmean(lats))).*(pi./180);
            xdistlon=xdist./(scalefactor);
            ydistlat=ydist./(111180);

          
            vel(isnan(dbzhc))=nan;
       
            figure(500);clf;set(gcf,'color','w','pos',[100 100 1500 500])
            subplot(1,2,1);cla;
            pcolor(xdist,ydist,vel);shading flat;% caxis([-10 20])
            caxis([-25 25])
            hold on;
            plot(0,0,'sk','markersize',10,'markerfacecolor','m')
            colormap(rbmapper(25,-25))
            xlim([-20 20]*10^3)
            ylim([-20 20]*10^3)



            box on; grid on;
            set(gca,'fontsize',15,'fontweight','bold','linewidth',2,'layer','top')
            daspect([1 1 1])
        end
    end
end

%% CONDUCT interpolation
Uint=FU(xdist(:),ydist(:),zdist(:)+alt(1));
Uint=reshape(Uint,size(xdist));
Vint=FV(xdist(:),ydist(:),zdist(:)+alt(1));
Vint=reshape(Vint,size(xdist));
%%
figure(100);clf; hold on;
subplot(1,3,1);cla; hold on;
pcolor(-xdist,ydist,Uint);shading flat;
caxis([-10 10])
colormap(rbmapper(10,-20))
subplot(1,3,2);cla; hold on;
pcolor(-xdist,ydist,Vint);shading flat;
caxis([-10 20])
colormap(rbmapper(10,-20))


%%
figure(500);clf;set(gcf,'color','w','pos',[100 100 1000 500])
hold on;
subplot(1,2,1);cla;hold on;

pcolor(xdist,ydist,vel);shading flat;% caxis([-10 20])
caxis([-25 25])
hold on;
plot(0,0,'sk','markersize',10,'markerfacecolor','m')
colormap(rbmapper(25,-25))
xlim([-5 20]*10^3)
ylim([-15 5]*10^3)
box on; grid on;
set(gca,'fontsize',15,'fontweight','bold','linewidth',2,'layer','top')
daspect([1 1 1])
xlabel('E-W Distance [m]');
ylabel('N-S Distance [m]')
cbh=colorbar;ylabel(cbh,'m/s')

figure(500);
subplot(1,2,2);cla; hold on;
vr=Uint.*sind(azmt).*cosd(evnow) + Vint.*cosd(azmt).*cosd(evnow);
pcolor(-xdist,ydist,-vr);shading flat;
caxis([-25 25])
hold on;
plot(0,0,'sk','markersize',10,'markerfacecolor','m')
colormap(rbmapper(25,-25))
xlim([-5 20]*10^3)
ylim([-15 5]*10^3)
box on; grid on;
set(gca,'fontsize',15,'fontweight','bold','linewidth',2,'layer','top');
daspect([1 1 1])
xlabel('E-W Distance [m]');
ylabel('N-S Distance [m]')
cbh=colorbar;ylabel(cbh,'m/s')


% export_fig DOW_v_WRF_RADIAL_WIND.png -m2.5


